/*
 * Assembler.cpp
 *
 *  Created on: 7 Sep 2011
 *      Author: Allan
 */

#include "Assembler.h"

// Boost includes
#include <boost/foreach.hpp>

// GeoReact includes
#include "Database.h"
#include "Phase.h"
#include "Multiphase.h"
#include "Reaction.h"
#include "ReactionSystem.h"
#include "System.h"
#include "Activity/AqueousActivity.h"
#include "Activity/GaseousActivity.h"
#include "Activity/MineralActivity.h"
#include "Kinetics/MineralKinetics/MineralKineticsDatabase.h"
#include "Utilities/Algorithms.h"
#include "Utilities/SpeciesUtils.h"

const System AssembleSystem(const vector<string>& species)
{
	vector<string> systemSpecies = Union(AssembleAqueousSpecies(species), species);
	
	vector<string> eSpecies = FilterSecondarySpecies(systemSpecies);
	
	Multiphase multiphase = AssembleMultiphase(systemSpecies);
	ReactionSystem eReactions = AssembleReactionSystem(eSpecies, multiphase);
	
	return System(multiphase, eReactions);
}

const System AssembleSystem(const vector<string>& species, const vector<string>& kSpecies)
{
	vector<string> systemSpecies = Union(Union(AssembleAqueousSpecies(species), species), kSpecies);
	
	vector<string> eSpecies = Difference(FilterSecondarySpecies(systemSpecies), kSpecies);
	
	Multiphase multiphase = AssembleMultiphase(systemSpecies);
	ReactionSystem eReactions = AssembleReactionSystem(eSpecies, multiphase);
	ReactionSystem kReactions = AssembleReactionSystem(kSpecies, multiphase);
	
	return System(multiphase, eReactions, kReactions);
}

const Multiphase AssembleMultiphase(const vector<string>& species)
{
	vector<Phase> phases;
	
	vector<string> aqueousSpecies = FilterSpecies(species, {"aqueous"});
	vector<string> gaseousSpecies = FilterSpecies(species, {"gaseous"});
	vector<string> mineralSpecies = FilterSpecies(species, {"mineral"});
	
	if(!aqueousSpecies.empty())
		phases.push_back(AssembleAqueousPhase(aqueousSpecies));
	
	if(!gaseousSpecies.empty())
		phases.push_back(AssembleGaseousPhase(gaseousSpecies));
	
	BOOST_FOREACH(const string& mineral, mineralSpecies)
		phases.push_back(AssembleMineralPhase({mineral}));
	
	return Multiphase(phases);
}

const Phase AssembleAqueousPhase(const vector<string>& species)
{
	return AssembleAqueousPhase(species, AqueousActivity(species));
}

const Phase AssembleAqueousPhase(const vector<string>& species, const AqueousActivity& activity)
{
	Phase phase;
	
	phase.SetName("Aqueous");
	phase.SetSpecies(species);
	phase.SetActivity(activity);
	
	return phase;
}
const Phase AssembleGaseousPhase(const vector<string>& species)
{
	return AssembleGaseousPhase(species, GaseousActivity(species));
}

const Phase AssembleGaseousPhase(const vector<string>& species, const GaseousActivity& activity)
{
	Phase phase;
	
	phase.SetName("Gaseous");
	phase.SetSpecies(species);
	phase.SetActivity(activity);
	
	return phase;
}

const Phase AssembleMineralPhase(const vector<string>& species)
{
	return AssembleMineralPhase(species, MineralActivity(species));
}

const Phase AssembleMineralPhase(const vector<string>& species, const MineralActivity& activity)
{
	Phase phase;
	
	phase.SetName("Mineral");
	phase.SetSpecies(species);
	phase.SetActivity(activity);
	
	return phase;
}

const ReactionSystem AssembleReactionSystem(const vector<string>& reactants, const Multiphase& multiphase)
{
	vector<Reaction> reactions;
	
	BOOST_FOREACH(const string& reactant, reactants)
		reactions.push_back(AssembleReaction(reactant, multiphase));
	
	return ReactionSystem(reactions);
}

const Reaction AuxiliarAssembleReaction(const ReactionInfo& rinfo, const Multiphase& multiphase)
{
	// The species that participate in the reaction
	vector<string> species = rinfo.products; 
	species.push_back(rinfo.reactant);
	
	// The stoichiometries of the species in the reaction
	vector<double> stoichiometries = rinfo.stoichiometries; 
	stoichiometries.push_back(-1.0);
	
	// Create the reaction equation
	ReactionEquation equation(species, stoichiometries);
	
	// Create the equilibrium constant of the reaction
	EquilibriumConstant K(rinfo.logKi, rinfo.Ti);
	
	// Create the reaction
	Reaction reaction(equation, K, multiphase);
	
	return reaction;
}

const Reaction AssembleReaction(const string& reactant, const Multiphase& multiphase)
{
	if(Database::GetInstance().GetSpeciesInfo(reactant).HasTag("basis"))
		return AssembleBasisReaction(reactant, multiphase);
	
	else if(Database::GetInstance().GetSpeciesInfo(reactant).HasTag("mineral"))
		return AssembleMineralReaction(reactant, multiphase);
	
	else 
		return AssemblePlainReaction(reactant, multiphase);
}

const Reaction AssemblePlainReaction(const string& reactant, const Multiphase& multiphase)
{
	// Get the information data of the reaction whose reactant is (reactant)
	const ReactionInfo reactionInfo = Database::GetInstance().GetReactionInfo(reactant);
	
	return AuxiliarAssembleReaction(reactionInfo, multiphase);
}

const Reaction AssembleBasisReaction(const string& reactant, const Multiphase& multiphase)
{
	ReactionEquation equation({reactant, reactant}, {-1.0, 1.0});
	
	Reaction reaction(equation, multiphase);
	
	return reaction;
}

const Reaction AssembleMineralReaction(const string& reactant, const Multiphase& multiphase)
{
	// Create the reaction object for the given mineral
	Reaction mineralReaction = AssemblePlainReaction(reactant, multiphase);
	
	// Look up at the database of mineral kinetics for the given mineral
	optional<MineralKineticsInfo> mineralKineticsInfo = 
		MineralKineticsDatabase::GetInstance().FindMineralKineticsInfo(reactant);
	
	// If the kinetics information of the minaral was found, create a mineral kinetics intance and set the mineral reaction with it
	if(mineralKineticsInfo)
	{
		// Create the mineral kinetics object of the given mineral
		MineralKinetics mineralKinetics(reactant, mineralKineticsInfo->mechanisms, mineralReaction, multiphase);
		
		mineralKinetics.SetSpecificSurfaceArea(mineralKineticsInfo->specificSurfaceArea);
		
		// Set the kinetics of the mineral reaction
		mineralReaction.SetKinetics(mineralKinetics);
	}
	
	return mineralReaction;
}

const vector<string> AssembleSpecies(const vector<string>& species)
{
	vector<string> primarySpecies = AssemblePrimarySpecies(species);
	
	vector<string> secondarySpecies = AssembleSecondarySpecies(species);
	
	return Union(Union(primarySpecies, secondarySpecies), species);
}

const vector<string> AssembleAqueousSpecies(const vector<string>& species)
{
	return FilterSpecies(AssembleSpecies(species), {"aqueous"});
}

const vector<string> AssembleGaseousSpecies(const vector<string>& species)
{
	return FilterSpecies(AssembleSpecies(species), {"gaseous"});
}

const vector<string> AssembleMineralSpecies(const vector<string>& species)
{
	return FilterSpecies(AssembleSpecies(species), {"mineral"});
}

const vector<string> AssemblePrimarySpecies(const vector<string>& species)
{
	// The set containing only unique basis species
	set<string> primarySpecies;
	
	// Make sure the hydrogen ion H[+] is present in the set of primary species in case H2O is listed in (species)
	if(IndexOf("H2O", species) < species.size())
		primarySpecies.insert("H[+]");
	
	// Loop over all the provided (species)
	BOOST_FOREACH(const string& s, species)
	{
		// Check if the current known species is a basis species in the database
		if(Database::GetInstance().GetSpeciesInfo(s).HasTag("basis"))
			primarySpecies.insert(s);
		else
		{
			// Get the info data of the reaction whose reactant is the current known species
			const ReactionInfo reactionInfo = Database::GetInstance().GetReactionInfo(s);
		
			// Add all the products in the reaction to the set of basis species being created (all these products are basis species in the database)
			primarySpecies.insert(reactionInfo.products.begin(), reactionInfo.products.end());
		}
	}
	
	// Convert the set of basis species to a vector container and return it
	return vector<string>(primarySpecies.begin(), primarySpecies.end());
}

const vector<string> AssembleSecondarySpecies(const vector<string>& species)
{
	const vector<string> primarySpecies = AssemblePrimarySpecies(species);
	
	// Get a reference to the list of reactions in the database
	const Database::ReactionInfoList& reactionInfoList = Database::GetInstance().GetReactionInfoList();
	
	// The vector of secondary species that can be built from the set of basis species
	vector<string> secondarySpecies;
	
	// Loop over all the reactions in the database
	BOOST_FOREACH(const ReactionInfo& r, reactionInfoList)
	{
		// An auxiliary boolean variable to determine if all the products of the current reaction belong to the set of basis species
		bool belong = true;
		
		//  Check if the all the products of the current reaction belong to the set of basis species (basisSpecies)
		BOOST_FOREACH(const string& product, r.products)
			if(!IsContained(product, primarySpecies))
				belong = false;
		
		// Insert the reactant (a secondary species) of the current reaction into the vector of secondary species
		if(belong)
			secondarySpecies.push_back(r.reactant);
	}
	
	return secondarySpecies;
}
